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PSEUDO  WIGNER-VILLE  DISTRIBUTION, 
COMPirrER  PROGRAM  AND  ITS  APPLICATIONS 
TO  TIME-FREQUENCY  DOMAIN  PROBLEMS 

by 

Jae-Jin  Jeon  and  Young  S.  Slun 


ABSTRCT 

Machinery  operating  in  non-stationary  mode  generates  a 
signature  wluch  at  each  instant  of  time  has  a  distinct  frequency. 
A  time-frequency  domain  representation  is  needed  to 
characterize  such  signature.  Pseudo  Wigner-Ville  distribution  is 
ideally  suited  for  portraying  non-stationary  signl  in  the  time- 
frequency  domain  and  carried  out  by  adapting  the  fast  Fourier 
transform  algorithm.  The  important  parameters  affecting  the 
pseudo  Wigner-Ville  distribution  are  discussed  and  sensitnity 
analyses  are  also  performed.  Practical  examples  of  an  actual 
transient  signal  are  used  to  illustrate  its  dynamic  features 
jointly  in  time  and  frequency. 


I.  INTRODUCTION 

The  physical  condition  or  state  of  health  of  machineries  which 
operate  in  transient  or  other  non-stationary  modes  are  difficult  to  predict 
with  any  degree  of  accuracy.  It  is  common  to  practice  periodic  preventive 
maintenance  on  ihese  machineries  in  order  to  avoid  failures  and  prolong  the 
useful  operating  life  of  the  equipment. 
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In  order  to  assess  the  physical  condition  of  machinery  without 
complete  disassembly,  a  physical  measurement  of  its  vibrations  is  conducted 
using  an  accelerometer.  Other  sensors,  such  as  temperature  or  pressure 
transducers,  could  also  be  used.  There  are  other  methods,  including  motor 
current  signature  analysis  on  electrically  driven  macliinery  and  wear  debris 
analysis  which  could  be  used.  However,  vibrations  are  used  predominantly 
for  machinery  condition  monitoring.  The  vibrations  are  recorded  in  the  time 
domain. 


There  is  a  need  for  a  method  to  represent  the  time  dependent  events 
which  occur  Nsdth  machinery  operating  in  non-stationary  modes.  At  each 
instant  in  time  as  the  speed  of  the  macliinery  changes,  the  frequency  content 
will  also  change.  The  pseudo  Wigner-V'ille  distributiontPWAT))  is  the  method 
which  was  chosen  to  portray  these  time  dependent  changes.  This  is  a 
continuation  of  work  initially  performed  and  published  by  Rossano,  Hamilton 
and  Shin  [Ij. 

The  pseudo  Wigner-Ville  distribution  is  a  three  dimensional  (time, 
frequency,  amplitude)  representation  of  an  input  signal  and  is  ideally  suited 
for  describing  transient  or  other  non-stationary  phenomena.  The  Wigner 
Distribution  (WDF)  has  been  used  in  the  areas  of  optics  [2,3,4]  and  speech 
analysis  [5.6].  Wahl  and  Bolton  [7]  used  it  to  identify  structure-borne  noise 
components.  Flandrin  et.  al.  [8]  recently  proposed  its  use  in  the  area  of 
machinery  condition  monitoring  and  diagnostics,  while  Forrester  [9]  is 
investigating  its  use  in  gear  fault  detection. 

For  such  a  non-stationary  signal  analysis,  spectrogram  is  commonly 
used,  which  is  based  on  the  assumption  that  it  is  a  collection  of  a  short 
duration  stationary  signals.  A  major  drawback  of  this  approach  is  that  the 
frequency  resolution  is  directly  affected  by  the  duration  of  short  stationary 
time,  wliich  subsequently  determines  the  time  resolution.  A  method  for  time- 
frequency  domain  signal  characterization  that  overcomes  this  drawback  is 
the  Wigner  distribution  which  was  first  introduced  by  Wigner  [10]  in  1932  to 
study  the  problem  of  statistical  equilibrium  in  quantum  mechanics.  The 


frequency  and  time  resolutions  of  the  Wigner  distribution  are  not  determined 
by  the  short  duration  but  rather  determined  by  the  selection  of  desired 
resolution  of  the  signal  itself. 

This  paper  discusses  the  important  parameters  affecting  the  P\V\'D 
in  order  to  machinery  condition  monitoring  and  presents  numerical  examples 
of  P\V'\'D  using  synthetically  generated  signals.  It  is  found  that  the  PW'VD 
is  very  effective  in  machinery  condition  monitoring  which  operates  in  non¬ 
stationary  modes. 


II,  PSEUnO  WIGNER-VILLE  DISTRII.  /IION 


A,  Wigner  Distribution  Function 

Signal  associated  with  most  Wbrational  phenomena  are  in  general  time 
varjdng,  which  means  that  their  characteristics  change  with  time  and  they 
have  various  features  in  difl'erent  time  frames.  The  general  spectrogram 
usually  requires  a  large  time-bandwidth  product  to  reduce  the  estimated  bias 
and  variability.  In  the  case  of  a  signal  containing  some  transients  or 
nonstationary  conditions,  the  traditional  approach  in  signal  analysis  fails  to 
describe  the  d>mamics  of  the  signal's  frequency  component  changes. 

The  general  expression  of  the  time-frequency  distribution  of  a  signal, 
w(t,(0)  is  given  by,  [11] 

\v(t,C0)  =  Ml  s''(u  - s(U  + dll  dl  d6  (1) 

2k  2  2 

where  s(u)  is  the  time  signal,  s*(.u)  is  its  complex  conjugate,  and  (P{0,  T]  is  an 

arbitrary  function  called  the  kernel.  By  choosing  different  kernels,  different 
distributions  are  obtained.  Wigner  distribution  is  obtained  by  taking  <p{6,  T) 

=  1.  The  range  of  all  integrations  is  from  -  co  to  oo  unless  otherwise  noted. 
Substituting  the  kernel  =  ItoEq.  (H, 

\v(t,a))  =  J■|^s''(u  --)  - 1)  s(u  +  -)  dx  du  (2) 

From  Eq.  (2l  the  Wigner  distribution  is  obtained, 

\v(t,to)  =  Js*(t  -^)  stt  -i-^)  dl  (3) 
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One  of  the  basic  frequency  representations  of  a  signal  is  the  power 
density  spectrum,  which  characterizes  the  signal's  frequency  component 
distribution.  The  power  spectral  density  function  S(CO)  of  a  signal  s(t)  can  be 
related  to  the  Fourier  transform  of  the  signal's  autocorrelation  function  R{T): 


Si(o)= 

with 

(4) 

A’(  r)  =  j  s(t)  s(t  +  r)di  . 

(51 

From  this  relation  a  time-dependent  power 
can  be  written  as 

spectral  density  function 

sv(i.{0)  =  J  RjiT) 

(G) 

where  now  R.(T)  is  a  time-dependent  or  local  autocorrelation  function,  Mark 
[12]  argued  for  symmetry, 

..  .  .  r  r. 

A\(r)  =  s  (i--)  s(i  +  -:^) 

which  gives  the  Wigner  distribution  function. 

B.  Properties  of  the  Wigner  distribution  function 

(7) 

The  properties  of  the  WDF  [13,14]  are  summarized  and  reinterpreted 
with  this  new  formulation  as  follows: 

(,i  )  The  WDF  is  a  real-valued  function. 

\v*(l,to)  =  wd.tj) 

(8) 
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(ii)  Tile  iiite^jral  of  the  WDF  with  respect  to  frequency  and  time  jields 
tile  instantaneous  signal  power  and  tlie  signal  s  power  density  spectrum, 
respectively. 

j wtt.to)  dm  =  2;r|s(t)i  • 

—  fo 

(9) 

|\v(t.m)dt  =  2;rjS(m)|  . 


(iii)  A  time  or  frequency  shift  in  the  signal  have  the  shift  in  the  WDF. 

Ifs(t)  ->  s(t  +  ),  then  W(t.m)  w(t  +  to.m),  tlO) 

ifs(l)  ‘s(t).  then  \v(t.m)  (0^)).  til) 

(iv)  The  WDF  is  symmetrica!  in  time  for  a  given  signal. 

Ifs(l)  — ^  st-t),  then\\tt,(U)  w(-t.m).  (12) 

ifs(()  s*(l).  then  w  tt.m)  ->  \v(t.-(0).  (13) 


(v  I  The  WDF  is  not  always  positive. 

(vi)  The  intogi'ation  of  the  square  of  the  WDF  equals  the  square  of  the 
time  integration  of  the  signal  s  power.  This  is  the  counterpart  of  l^arseval's 
relation  of  the  WDF,  called  Mayol's  fomula. 


J  |\v(t.m);  dt  dm  = 


J  s“(l )  dl 


(14) 


(vii)  The  WDF  possesses  basic  symmetry  under  the  interchange  of  time 
and  frequency  parameters  with  the  Fourier  transform  of  a  given  signal. 

Let 

1 

s(i)  =  —  I  c  .S(  m)  dm;  (15) 

2/:  ’ 
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then 


w(t,w)  =  “ 

2n 


j  S(4)  +  S*{(0  -  ~)  d^, 


(16) 
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III.  IMPLEMENTATION  WITH  DIGITAL  SIGNAL  PROCESSING 


A.  Discrete  Wigner  Distribution  Function 

There  are  two  distinct  advantages  for  the  calculation  of  the  WDF. 
First,  it  has  the  form  of  the  Fourier  transform  and  the  existing  FFT 
algorithm  can  be  adapted  for  its  computation.  Second,  for  a  finite  time  signal, 
its  integration  is  finite  wnthin  the  record  length  of  the  existing  signal. 

The  discrete  time  Wigner  distribution  as  developed  by  Claasen  and 
Mecklenbrauker  [13]  is  expressed  by, 

w'(t,co)  =  2  X  +  x)  s'"!!  -  x)  (17) 

The  discrete  version  of  Eq.  (17)  for  a  sampled  signal  s(n),  where  n=0  to  N-1, 
has  the  form, 

1  ^'“1  * 

w(7,k)  =  —  s{(  +  n)  s  {( -n)  c  '  ,  k=0,l,2,...N-l  (18) 

Nn=0 


where  s(m)=0  for  m  <  0  and  ni  >  N-1.  However,  in  order  to  utilize  the  FFT 
algorithm,  it  must  be  assumed  that  the  local  autocorrelation  function  has  a 
periodicity  of  N.  This  is  just  for  operational  convenience  and  should  not  apply 
to  the  interpretation  of  s(m)  Eq.  (18)  can  be  rewritten  as, 

1  N'l  *  -y^n(k  +  m— ) 

w[f,k -h  m(N  /  2)]  =  —  +  (/-n)e  ^ 

^  n=0 

1  N-1  -i— uk 

=  —  y  s(7-»-n)  s*((-n)  e  (19) 

N 

n=0 

=  w(Gk) 
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since  =  1  for  m=integers. 

Eq.  (19)  indicates  that  the  WDF  has  a  periodicity  of  N/2.  Hence,  even 
when  the  sampling  of  s(t)  satisfies  the  Nyquist  criteria,  there  are  still  aliasing 
components  in  the  WDF.  A  simple  approach  to  avoid  aliasing  is  to  use  an 
analytic  signal  before  computing  the  WDF.  In  1948,  J.  Ville  [15]  proposed  the 
use  of  the  analjlic  signal  in  time-frequency  representations  of  a  real  signal. 


B.  Analjdic  Signal 

An  analytic  signal  is  a  complex  signal  which  contains  both  i '  al  and 
imaginary’  components.  The  advantage  of  using  the  analytic  signal  is  that  in 
the  frequency  domain  the  amplitude  of  negative  frequency  components  are 
zero.  This  satisfies  mathematical  completeness  of  the  problem  by  accounting 
for  all  frequencies,  yet  does  not  limit  the  practical  application  since  only 
positive  frequency  components  have  a  practical  interpretation.  The 
imaginary  part  is  obtained  by  Hilbert  transform.  The  analytic  signal  may  be 
expressed  by, 

S(t)  =  Sr(t)  +  j  H{Sr(t)}  (20) 

where  H(Sr(t)}  is  a  Hilbert  transform  and  generated  by  the  convolution  of 
the  impulse  response  h(t)  of  a  £0-degree  phase  shift  as  follows: 

H{Sr(t)}  =  Sr(t)  *  h(t) 

h(t)  =  2  sin-(;n/2) 

m 

=  0, 

where  *  denotes  the  convolution.  Rewriting  Eq.  (21)  to  discrete  form, 

oo 

H{Sr(n)}  =  Xh(n -m)  Sr(m)  (22) 

m=-oo 


(21) 

t  =  0, 
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The  distribution  resulting  from  an  analytic  signal  being  processed  through 
the  Wigner  distribution  is  commonly  termed  as  Wigner-Ville  distribution. 


C.  WDF  with  Digtal  Signal  Processing 


To  calculate  the  Wigner  distribution  of  the  sampled  data,  it  is 
necessary  that  Eq.  (18)  be  modified  to  Eq.  (23),  because  the  WDF  has  N/2 
periodicity. 


2N-1 


w(rnAt,kA(y)  =  2At  J^s[(ni  +  n)At]  s* [(m  - n)At]  ^ 

n=0 


-y2/nik/(2N)  (23) 


where  A(0  =  71  /  (2NAt)  and  At  is  the  sampling  interval.  The  algorithm  used 
in  this  paper  is  based  on  one  written  by  Wahl  and  Bolton[7]  and  can  be 
expressed  as: 

w(niAi,kAa))  =  Re  [2At  FFT(corr(i))] 

corT(i)  =  s(m+ i-1)  s*(m-i  +  l),  m  >  i  (24) 

=  0,  m  <  i 

where  1  <  i  <  N  +  1 

corr(2N  -  i  +  2)  =  corr  (i),  2<i<N 


The  frequency  resolution,  A(0 ,  in  Eq.  (23)  is  different  from  that  obtained  by 
FFT  of  the  original  N  point  time  record  in  two  respects.  The  first  difference  is 
that  the  argument  of  the  time  signal  and  its  conjugate  contains  a  factor  of 
1/2,  and  secondly,  the  autocorrelation  of  the  time  signal  is  twice  the  length  of 
the  original  record  and  therfore  the  FFT  is  evaluated  over  2N  points.  The 
result  is,  that  the  WDF  frequency  resolution  is  one  forth  the  resolution  of  an 
ordinary  power  spectrum  density  function. 


Before  processing  the  WDF,  a  modified  Hamming  window  is  applied  to 
the  time  domain  signal  to  reduce  the  leakage  caused  by  the  discontinuity  of 
the  finite  record  of  data,  which  will  be  called  as  data  tapering.  This  type  of 
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,/ 

window  is  preferable  since  it  alters  the  amplitude  of  fewer  data  points  at  the 
beginning  and  the  end  of  the  data  block.  A  modified  Hamming  window,  D(t)  is 
given  by: 

0.54  -  0.46  *  cos(  10;:t/r),  0  <  t  <  T/10, 

D(t)  =  {  1.0,  T/10  <  t  <  9T/10,  (25) 

0.54  -  0.46  *  cos(10rt(T-t)/T),  9T/10  <  t  <  T. 


D(t) 


t - ► 

Fig.  1  Modified  Hamming  Window 


Two  other  characteristics  of  the  WDF  should  be  also  noted.  First,  the 
WDF  of  the  sum  of  two  signals  is  equal  to  the  sum  of  the  WDF  of  each  signal 
plus  cross  term  that  appear  when  the  cross-correlation  of  the  two  signal  is 
non-zero.  Second,  the  WDF  may  have  negative  values,  which  may  be  largely 
caused  by  interference  due  to  the  presence  of  these  cross  terms.  In  the  case 
of  input  signals  that  contain  multi-frequency  components,  the  Wigner-Ville 
distribution  of  most  signals  are  very  complicated  and  difficult  to  interpret. 

There  are  two  methods  to  suppress  the  interference  components  of  the 
WDF.  Claasen  and  Mecklenbrauker[12]  describe  the  application  of  a  sliding 
window  in  the  time  domain  before  calculating  WDF.  The  WDF  obtained  with 
a  window  function  is  called  the  Pseudo-  Wigner  distribution  function.  A 
second  option  is  to  smooth  the  WDF  with  a  sliding  averaging  window  in  time- 
frequency  plane.  In  both  case  the  result  is  to  deemphasize  components 
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arising  from  calculations  and  to  emphasize  deterministic  components. 
Obviously,  averaging  a  Wigner-Ville  distribution  will  result  in  a  Pseudo 
Wigner-Ville  distribution. 


In  this  research,  a  sliding  exponential  window  in  the  time-frequency 
domain  was  chosen.  That  is,  a  Gaussian  window  function,  G(t,  co)  is  selected 
to  reduce  the  interference  and  to  avoid  the  negative  values  as  follows: 
let 


G(t.ca) 


then 


1 


co“  ^ 

-( - T  +  — r) 


la, 


ia 


27rOtOto 


(0 


(26) 


\v(t.(0) 


—  G(t  -  t'.co-to')  di'dco'  >0 

2k" 


(27) 


where  Gp  >  0  and  >  1/2  [16].  The  time  and  the  frequency 

resolution  At  amd  Aco  of  this  Gaussian  window  are  related  by. 

Of  =  j  At,  G<.^  =  k  Aw  (28) 

in  the  discrete  form.  Then  the  condition  for  the  WDF  to  be  positive  in  this 
case  is 

j  Ai  k  ACO  >  1/2.  (29) 

This  is  the  time-frequency  version  of  Heisenberg's  uncertainty  relation[14]. 
If  the  segmentation  of  time  and  frequency  for  a  given  signal  from  Eq.  (3) 
violates  this  uncertainty  principle,  the  corresponding  WDF  may  not  be 
positive. 

To  perform  the  convolution  on  the  sampled  \\T)F,  the  Gaussian  window 
function  was  applied  to  the  range  ±2G{  and  ±2  0^,.  Selecting  co  and  t  to  be 

the  multiple  of  time  and  frequency  steps,  the  sampled  Gaussian  vcdndow' 
function  is  expressed  by. 
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(30) 


G(p.q) 


271  j  k  AlAco 


where  p  and  q  are  an  integer  numbers  in  the  range  ±2j  and  ±2k,  respectively. 
The  convolution  of  the  sampled  \\T)F  and  the  Gaussian  wdndow  function  can 
be  evaluated  as  follows: 

At  A(0 

v.''(7',m)  =  — -  1  Zv.(p.q)  G(p- q -m)  (31) 

p=/'-j  q=m-k 

whore  \v'(/,m)  is  the  smoothed  \VT)F  or  Pseudo  Wigner-Ville  distribution. 


Fig.  2  shows  a  block  diagram  for  computational  sequence  of  the  Pseudo 
Wigner-Ville  distribution.  A  time-varjdng  signal  sampled  with  the  Nyquist 
rate  is  first  high  passed  through  a  digital  filter  if  the  signal  involves  the  zero 
frequency  component,  i.e.,  DC  component,  and  converted  into  the  analjlic 
signal  through  a  Hilbert  transform.  Then,  the  time-dependent  correlation 
function  is  computed  and  the  result  is  the  WTDF  in  terms  of  both  time  and 
frequency  domain  by  FFT.  The  final  step  is  to  i\>mpute  the  convolution  with  a 
Gaussian  v-indow. 


Fig.  2  Computational  block  diagram  of  Pseudo  Wigner-Ville  Distribution 
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D.  Highpass  Digital  Filter 


Filters  are  a  paticularly  important  class  of  linear  time-invariant 
system.  Strictly  speaking,  the  term  frequency-selective  filter  suggests  a 
system  that  passes  certain  frequency  components  and  totally  rejects  all 
others,  but  in  a  broader  context  any  system  that  modifies  cetain  frequencies 
relative  to  others  is  also  call  filter.  The  design  of  filters  involves  the  following 
stages  :  (1)  the  specification  of  the  desired  properties  of  the  system;  (2)  the 
approximation  of  the  specifications  using  a  causal  discrete-time  system;  and 
(3)  the  realization  of  the  system. 


In  this  paper,  Nonrecursive(finite  impulse  response  -  FIR)  highpass 
filter  was  used  for  the  elimination  of  undesired  low  frequency  components. 
The  basic  design  was  to  use  a  symmetric  filter  of  the  form, 


Nt 


y(i)  =  X  ^k  s(‘ 

(32) 

k  =  -M 

with 

b-k  =  bk 

(33) 

and 

^  sin  2TBkT 

bk  = 

(34) 

;rk 


where  b^.  is  the  filter  w'eights,  yti)  is  the  filtered  signal,  B  is  the  cutoff 
frequency,  T  is  a  sampling  interval  and  M  is  the  span  of  the  filter;  2M+1 
weights  are  employed  because  of  ssTninetry.  only  M-fl  need  be  generated. 

The  b|^  weights  are  computed  over  the  range  -M  to  M.  The  weights  are 

multiplied  by  a  window  function.  Potter  discusses  a  number  of  window's  in 
the  referenced  work.  His  P310  window  was  found  to  be  appropriate  for  filter 
implementation.  It  takes  the  form. 


w 


k 


w'here 


-t 


0 


3 


I  dp 

p=.3 


COS 


;rpk 


(35) 
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and 


k  =  ±  M 
Otherwise 


do  =  1 

d.i  =  d,  =  0.684988 
d.2  =  d.  =  0.202701 
d.3  =  d,  =  0.0177127 


vv  =  do  +  2  X  d  =  2.8108034 

P-3 


(36) 


^37) 


For  a  highpass  filter  with  pass  band  from  the  cutof['frequency(B)  to  the 
maximum  frequency,  generate  a  low  pass  filter  on  the  range  C  -  B,  and  then 
subtract  the  central  weight  from  unity  and  change  the  signs  of  the  remainder 
of  the  weights. 
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IV.  EXAMPLES  AND  DISCUSSIONS 


Maclunery  operating  in  transient  mode  generates  a  signature  in  which 
the  frequency  content  varies  at  each  instant  of  time.  To  characterize  such 
signatures  and  to  understand  the  vabrational  behavior  of  such  machineries, 
time-frequency  domain  representation  of  the  signal  is  needed.  As  discussed 
in  the  previous  sections,  Wigner  distribution  is  a  signal  transformation  that 
is  particularly  suited  for  the  time-frequency  analysis  of  nonstationary’  signals. 
There  are  many  advantages  of  using  P\WD  for  both  steady  and  transient 
signals.  However,  there  are  also  several  disadvantages,  for  example,  the 
drastic  increase  of  peak  value  when  the  frequency  content  of  signal  changes 
abruptly.  A  computer  program  has  been  developed  for  PWVD.  Two  different 
versions  are  available  at  the  present  time;  workstation  and  IBM  PC 
compatible. 


A.  Harmonic  Wave 

Fig.  3  shows  the  PWTD  of  the  pure  sine  wave  with  two  frequency 
components  (lOOHz,  400Hz),  respectively.  The  modified  Hamming  window 
was  applied  to  the  time  domain  signal  and  the  Gaussian  smootlhng  window 
function  was  applied  on  time-frequency  domain  Winger-Ville  distribution. 
The  slope  of  the  end  edges  are  due  to  data  tapering  by  using  the  modified 
Hamming  window.  Fig.  4  shows  the  PWAD  of  the  sine  wave  that  have  the  10 
and  50  signal  to  noise  ratio,  respectively.  The  shape  of  P\V\^D  is 
changed  at  the  crest  by  the  contamination  of  noise.  The  crest  has  the 
complicated  shape  with  decreasing  of  signal  to  noise  ratio.  However,  the 
PWVD  well  represents  the  signal  components  from  the  given  signal  with 
noise.  The  practical  example  is  shown  in  Fig.  14.  The  notation  fg  and  N 
used  in  the  figures  are  sampling  frequency  and  the  total  number  of  time  data 
points. 
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'o  6.0 
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Fig.  3  Pseudo  Wiguer-Ville  distribution  of  100  and  400  Hz  Pure  Sine  Waves 
(fg=1000  Hz,  N=256  and  Smoothing  Window  Size=10xl0) 
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Fig.  4  Pseudo  Wigner-Ville  distribution  of  300  and  750  Hz  sine  waves;  signal 
to  noise  ratio  (a)  10  9c  and  (b)  50  9r  . 

(fs  =  2048  Hz,  N=1024  and  smoothing  windosv  size=18xl8) 
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B.  Harmonic  Wave  with  Stepwise  Frequency  Changes 

Fig.  5  shows  the  P\SVD  of  the  sine  wave  with  500  Hz  in  the  time  from 
0.085  sec  to  0.17  sec.  The  PVV'VT)  well  represents  the  time  delay  of  the  signal. 
The  Fig.  6  shows  (a)  the  sine  wave  with  stepwise  frequency  changes,  100  Hz, 
250  Hz  and  500  Hz  and  (b)  its  P\V'\T).  The  P\S'\T)  shows  the  time  delay  and 
frequency  component  of  the  signal.  The  wide  spread  of  PWTD  at  the  edge  of 
each  frequency  region  is  noticed.  This  phenomenon  is  caused  by  the 
discontinuity  of  the  signal  in  time  domain  and  the  leakage  in  digital  signal 
processing.  This  effect  may  be  reduced  by  applj-ing  the  data  tapering  to  the 
actual  signal  block.  Nevertheless  the  PWWD  represented  the  characteristics 
of  the  signal  well.  PW\’D  can  portray  the  charactenstics  of  the  steady  state 
signals  involving  time  delay  and  multi-frequency  components  .  If  diflerent 
size  of  the  smoothing  window-  are  applied,  the  PWVD  amplitude  changes,  but 
the  total  energy  remains  unchanged. 

0.0^00 


0.007® 


Fig.  5  500  Hz  sine  wave  with  finite  duration 
(fg=2000  Hz,  N=512  and  Smoothing  Window  Size=10xl0) 
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C.  Composite  Signal  with  Two  Frequency  Components  at  Each  Time 

The  P\S'\T)s  of  the  nons.  ationary  signals  were  studied  and  the  results 
were  shown  in  Fig.  7  through  10.  Fig.  7  shows  (a)  the  time  signal  composed 
of  two  sweeping  frequency  components  at  each  time,  one  increasing  and  the 
other  decreasing  with  the  same  rate,  and  (b)  its  Wigner-Ville  distribution 
(before  appljnng  the  smoothing  ssindow)  (c)  its  contour  plot  of  WDF  and  (d)  its 
pseudo  Wigner-Ville  distribution  (after  applying  the  smoothing  window), 
respectively. 


The  efTect  of  cross  (or  interference)  term  is  significant  and  appeared  in 
the  average  frequency  region.  This  is  one  of  the  disadvantages  of  using 
Wigner-Ville  distribution  but  it  is  a  characteristic  of  the  distribution.  W’hen 
Gaussian  window  was  applied  to  Wigner-Ville  distribution,  the  effect  of  cross 
term  disappeared.  The  main  lobe  of  PWVD  is  wider  and  its  amplitude  is 
significantly  reduced.  The  large  peak  at  the  intersection  point  of  two 
sweeping  frequency  signals  is  mainly  caused  by  the  doubling  effect  of 
amplitudes  of  two  signals. 


c 

X) 

a 

♦.» 

E 

09 

f) 


10.0 


s.o 


0.0 


-6.0 


•10.0 


0.0  0.2  0.4  O.B  0.6  1.0 


Time  (sec) 
(a)  Time  signal 
Fig.7  (continued) 
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(d)PWVD 


Fig.  7  Composite  Signal  with  Two  Frequency  Components  at  Each  Time 

s(t)=4cos(2Tc  32t^)  +  4  cos{2rt(40+32(2-t)]t} 

(f  =256  Hz,  N=256  and  Smoothing  Window  Size=10xl0) 


D.  Linear  Chirp  Signal 


Another  type  of  a  non-stationary  signal  sweeps  up  and  down  in 
frequency  is  called  a  linear  chirp  signal  and  is  shown  in  Fig.  8(a).  This  signal 
has  only  one  frequency  component  at  each  time.  The  effect  of  cross  terms 
appears  in  the  Wigner-Ville  distribution,  as  shown  in  Fig.  7(b).  The 
smoothing  window  was  applied  to  Wigner-Ville  distribution  and  the  result  is 
shown  in  Fig.  7(d).  Fig.  8  (c)  is  the  contour  plot  of  WDF.  As  expected,  the 
effect  of  cross  term  is  significantly  reduced.  However,  the  unusual  peak 
(called  'ghost'  peak)  appeared  at  the  point  where  the  direction  of  sweep 
changes.  To  understand  the  cause  of  this  phenomenon,  the  PWW^D  was 
integrated  along  the  frequency  axis  and  it  was  found  that  the  square  root  of 
the  resultant  amplitude  was  the  amplitude  of  original  time  signal,  implying 
that  the  energ>’  content  remained  constant.  The  following  function  was  used 
to  generate  the  linear  chirp  signal; 


s(t)=  sin 
s(t)=  -sin 


2;r  30  + 


2m 


30 


220(i-l)3 
256  J  \ 

220(512-i)^ 
256 


(0.256 -t) 


1  <  i  <  256 

256  <  i  <  512 


(38) 


where  t  =  (i-1)  dt  and  dt=0.0005. 


Fig.8  (continued) 
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(c)  Contour  plot  of  WDF 
Fig.8  (continued) 


E.  Composite  Signal  of  Sweeping-up  and  Steady  Frequency 

The  signal  which  sweeps  up  along  the  frequency  for  first  0.5  second 
and  holds  to  a  constant  frequency  for  next  0.5  second  was  considered.  This 
signal  is  typical  speed  profile  of  start-up  stage  of  pump.  Fig.  9  shows  (a) 
PWVD  and  (b)  its  contour  plot.  The  interesting  phenomenon  was  observed  in 
PWVD  that  the  sweep-up  portion  of  signal  (first  half  seconds)  has  a  lower 
amplitude  and  wider  main  lobe  compared  with  the  steady  fi-equency  region  of 
signal  (second  half  seconds).  When  the  PWVD  was  integrated  along  the 
frequency  axis  and  it  was  found  that  the  resultant  amplitudes  in  these  two 
regions  are  same.  The  following  functions  were  used  to  generate  the  desired 
signal: 


s(t)  =  4cos(27t32t“ ),  0<t<0.5sec. 

s(t)  =  4cos(27t64t),  0.5  <  t  <  1.0  sec. 

Fig.  10  is  the  PWVD  of  the  signal  which  sweeps  up  along  the  frequency 
with  a  logarithmic  rate,  that  is,  the  sweep  rate  is  propotional  to  the  square 
root  of  time.  It  was  found  that  the  maximum  magnitude  of  the  PWVD 
increases  with  increasing  the  frequency.  This  fact  is  shown  that  the  PWVD 
of  the  stable  signal  has  a  larger  magnitude  than  the  unstable  signals 
although  having  the  same  magnitude  in  time  domain  and  the  PWVD  is  the 
good  tool  for  the  analysis  of  the  stability  of  the  signal.  The  following 
functions  were  used  to  generate  the  desired  signal: 

s(t)  =  4cos{2;r(30  +  60  t''^‘)t}.  (40) 

The  instantaneous  frequency  of  the  Eq.(40)  is  the  derivative  of  the  argument 
of  the  cosine  function  which  is 

/(t)  =  90  -f-  30  (41) 
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Tim*  (s*c) 


(b)  Contour  plot  ofPVWD 

Fig.  9  of  a  Composite  Signal  of  Sweeping-up  and  Steady  Frequency 

(f  =256  Hz,  N=256  and  Smoothing  Window  Size=10xi0) 


F.  Sweep  Rate  Effect 


The  effect  of  sweep  rate  on  PW'N’D  was  investigated.  The  sweep  rate  is 
the  frequency  change  per  unit  time.  The  power  spectrum  density  of  a  U^pical 
swept  sine  is  shown  in  Fig.  11.  The  incorrect  assumption  is  often  made  that  a 
swept  sine  of  constant  amplitude  has  a  flat  spectrum.  As  can  be  seen  from 
the  plot  this  is  not  so.  Fig.  12  shows  the  PW'VTis  of  the  linear  chirp  signal 
with  a  various  sweep  rates:(a)  has  zero  sweep  rate  and  (b)  has  lower  sweep 
rate  than  (c).  It  can  be  seen  that  the  amplitude  of  PWTD  decreases  with 
increasing  sweep  rate  but  energ>’  remains  unchanged.  This  result  appeared  to 
be  caused  by  Heisenberg's  uiicertainty  relation  between  time  and  frequency. 
However,  based  on  this  study,  it  is  clear  that  the  'ghost'  peak  (see  Fig.  8) 
appears  due  to  the  instantaneous  zero  sweep  rate  at  the  point  where  the 
direction  of  ssveep  changes.  Also  the  peak  value  is  affected  by  the  size  of 
smoothing  window. 
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Fig.  11  Swept  sine  wave  spectimm. 
(Frequency  range  1  to  20  Hz,  sweep  time  22.5  sec) 
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G.  Haniionic  Wave  with  Some  Glitches 

The  interesting  phenomena  on  the  signal  with  an  abnormal 
components  as  a  fault  were  investigated.  Fig  13  shows  the  PWVD  of  the 
harmonic  wave  with  glitch  at  a  small  region  of  the  time  record:  (a)  is  the  time 
signal,  (b)  is  the  PWAT)  and  (c)  is  the  contour  plot  of  the  PW\T).  It  can  be 
seen  that  the  PW^T)  of  the  signal  Fig.  13(a)  well  represents  the  loacation  of 
each  glitch  and  its  frequency  components.  This  characteristic  of  the  PWWD  is 
useful  to  detect  the  faults  or  glitch  and  to  monitor  the  condition  on  the 
vibrational  machinery  having  the  periodicity  such  as  a  gear  train.  The 
general  rotating  machinery  has  a  periodic  signal  pattern  in  time  domain. 


0.00  0.26  0.60  0.76  1.00 

Tlm«  (sec) 


Fig.  13  (continued) 


H.  Actual  t'an  Signal 


The  acceleration  signal  of  a  fan  was  measured  at  the  steady  state 
speed  and  the  result  was  shown  in  Fig.  14.  The  crest  has  the  complicated 
shape  on  time  axis  as  shown  in  Fig.  4.  The  first  peak  is  the  fundamental 
frequency  of  the  blade  rate  and  the  second  peak  is  3rd  harmonics.  The  third 
peak  is  the  fundamental  frequency  of  motor  by  the  pole.  The  measured 
vibration  signal  was  contaminated  with  the  noise.  If  the  measured  signal 
involves  the  faults,  the  PWVD  sviW  represent  the  different  pattern  having  the 
abnormal  frequency  components  in  comparison  with  the  normal  condition 
with  time. 


Fig.  14  PW'MD  of  the  actual  fan  signal  in  the  constant  speed, 
(fs  =  512  Hz,  N=512,  smoothing  window  size  13x13) 


I. 


Actual  Pump  Start-up  RPM  Signal 


The  start-up  transient  speed  of  the  pump  was  measured  and  the 
results  were  shown  in  Fig.  15.  The  time  signal  is  measured  by  magnetic 
sensor  and  has  the  same  magnitude  independently  on  time.  The  PWVD  is 
shown  in  Fig.  15(a)  and  the  contour  view  is  shown  in  Fig.  15(b).  The  contour 
plot  shows  that  the  speed  of  the  pump  runs  up  when  initially  started,  reaches 
the  maximum  RPM  and  coasts  down  gradually.  Near  the  maximum  speed 
during  the  run  up,  the  sweep  rate  was  rapidly  decreased  and,  as  a  result,  the 
peak  value  was  rapidly  increased.  When  the  sweep  rate  is  close  to  zero  at  the 
normalized  time  of  0.4,  the  amplitude  attains  the  maximum  value. 


Fig.  15  (continued; 
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Fig.  15.  Pseudo  Wigner-Ville  Distribution  of  Transient  Speed  of  the  Pump; 
(a)  PWVD  and  (b)  contour  plot. 
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V.  CONCLUSIONS 


The  pseudo  Wigner-Ville  distribution  has  been  investigated  and  applied 
to  analyzing  non-stationary  signals  t>’pical  of  transient  machinery'  signatures. 
The  results  of  this  research  will  be  a  valuable  asset  for  condition  monitoring 
of  transient  machinery-.  The  following  conclusions  can  be  drawn: 

(1)  The  pseudo  Wigner-Ville  distribution  is  ideally  suited  for  portrajdng  non- 
stationar}’  time  signals. 

(2)  The  use  of  modified  Hamming  window  to  time  signals  is  effective  to 
reduce  the  edge  effect  of  discontinuity. 

(3)  The  use  of  the  analytic  signal  in  calculating  the  Wigner  distribution 
eliminates  aliasing  problem. 

(4)  The  Gaussian  window  function  for  smoothing  the  Wigner-Ville 
distribution  is  very  effective  and  the  presence  of  cross  terms  is  significantly 
reduced. 

(5)  Both  the  amplitude  and  the  main  lobe  of  the  pseudo  Wfigner-Ville 
distribution  is  significantly  affected  by  the  sweep  rate.  As  the  absolute  sw'eep 
rate  increases,  the  amplitude  of  the  PWVD  decreases  and  the  main  lobe 
becomes  wider. 

(6)  The  PWWD  characterizes  the  time-frequency  domain  distribution  of  the 
signal  well  and  may  be  useful  tool  for  the  machinery'  condition  monitoring. 
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APPENDIX  A.  USER'S  GUIDE  OF  THE  PROGRAM 


This  program  calculates  either  the  Wigner-Ville  distribution  or  the 
smoothed  Wigner-Ville(Pseudo  Wigner-'Ville)  distribution(PWVD)  of  a  time 
series.  This  program  uses  the  FORTRAN  77  language  and  it  is  possible  to 
run  at  workstation  and  IBM-PC  compatible  computer.  The  user  supplies 
the  real  data  of  the  time  histor>’,  the  sampling  parameters,  the  digital  filter 
parameters,  the  smoothing  parameters  and  the  output  parameters.  The 
program  outputs  2-D  array  containing  the  Wiger-  'Ville  distribution 
(rwdf.out)  and  the  smoothed  Wigner-Ville  distribution,  i.e,  the  pseudo 
Wigner-Ville  distribution  (rswdf.out)  dependently  output  parameters. 


A-1.  PROGRAM  INPUT/OUTPUT  (Main  Program) 

This  program  is  interactive  and  displays  an  explanations  about  a 
desired  values.  Fig.  A-1  is  a  flowchart  of  this  program  and  the 
computational  block  diagram  is  shown  in  Fig.  2. 

1)  INPUT 

This  program  displays  the  input  variables  as  follow. 

Enter  name  of  signal  input  file 

The  user  must  supply  the  input  filename  under  25  characters. 

Number  of  the  sampled  data  point 

The  user  must  input  the  number  of  sample  data(dp).  If  dp  is  larger 
than  2048,  the  user  must  change  np  in  a  available  memory  size  of  the  user's 
computer  and  recompile  the  source  program.  The  dp  must  be  a  multiple  of 
2  because  of  FFT. 


Fig.  A-1  Flow  chart  (continued) 


44 


The  next  step  is, 


Do  you  wish  to  remove  the  mean  value  ? 

Enter  1  for  yes  or  0  for  no. 

The  user  inputs  the  desired  value  in  according  to  the  characteristics 
of  input  time  signal.  If  the  time  signal  have  the  DC  components,  the  user 
had  better  select  mvopt=] . 


Do  you  want  to  apply  highpass  digital 
filter  to  the  original  data  ?  (y/n) 


If  'y'  is  selected,  the  desired  cutoff  (half-power  point)  frequency  is  inputed 
after  displaWng  the  follow  command. 

Enter  the  cutoff  frequency  of 
the  digital  highpass  filter  (in  Hz) 

This  program  uses  the  non-recursiveffinite  impulse  response  -  FIR) 
highpass  filter  for  the  elimination  of  undesired  low  frequency  components. 
If  the  user  wants  to  pass  the  highpass  filter,  input  the  desired  cutoff 
frequency  in  Hz. 

The  next  procedure  is, 


Input  tlie  desired  reduction  size 
input  1  for  64  by  32 
input  2  for  128  by  64 
input  3  for  128  by  128 
input  4  for  256  by  128 
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This  input  statement  describes  the  desired  output  file  size  for  the  output 
and  the  plotting  because  the  original  VVigner  distribution  involving  the  all 
domains  needs  the  huge  memory  size. 

As  the  next  step,  tliis  program  chooses  the 

The  desired  smootfhng  window  size 

input  the  smoothing  parameter  for 
frequency  in  Gaussian  funclion(nO 

input  the  smootliing  parameter  for 
time  in  Gaussian  function(mt) 


The  smoothing  parameters  for  frequency  (nO  and  time  (mt)  depend  on  the 
number  of  the  sampled  data.  The  user  inputs  the  integer  values  satisfying 
the  following  condition. 


nf  *  nit  > 


n 


tA-i) 


We  recommend  to  select  the  square  smoothing  \s'indow  size. 


2)  OUTPUT 

This  program  generates  the  two  output  files.  One  is  the  file  for  the 
unsmoothed  Wigner-Ville  distribution  and  the  other  is  the  file  for  the 
smoothed  WignerA’ille  distribution.  The  output  files  consist  of  one  column 
along  with  the  time  and  frequency  axis  dependently  the  output  parameter 
as  followng  descriptions. 


O.OOOOOE+OO 

O.IOOOOE+Ol 


:  initial  time 
:  time  record  length 


O.OOOOOE+OO 

0.25600E+03 

0.53830E-06 

0.49456E-06 

0.18271E-0G 


:  starting  frequency 
:  end  frequency 
64  ;  reduction  size 


:  the  results  of  PWVD  or  W'VT) 


The  graphic  output  of  the  results  used  of  CA-DISSPLA  version  11.0. 
The  3-D  graphic  program  describes  the  unsmoothed  or  smoothed  Wigner- 
Ville  distributionior  pseudo  Wigner-Ville  distribution).  The  2-D  graphic 
program  describes  the  contour  plot  of  the  results  and  the  time  series. 


A-2.  EXPLANATIONS  OF  INPUT  VARIABLES 


inname  Input  file  containing  the  real  time  signal, 
dp  The  number  of  the  sampled  data  point, 

mt  Smoothing  parameter  for  time  in  Gaussian  Fn. 

nf  Smoothing  parameter  for  frequency  in  Gaussian 

Fn. 

mvopt  Option  with  respect  to  remoWng  mean  value 
if  mvopt  =  1,  then  zero  mean, 
if  mvopt  =  0,  then  no. 

bw  Cutoff  frequency  of  highpass  filter 


4  9 


A-3.  SUBROUTINES 


1)  indata(dp,tin,ain) 

This  subroutine  reads  the  input  file  'inname'.  The  input  data  file 
involves  the  time  and  amplitude  in  real  value.  The  reading  format  is  a  free 
format. 


2)  dtcalc(dp,tin,dt) 

This  subroutine  calculates  the  mean  time  interval  dt  from  the  input 
time  signal.  Because  the  obtained  data  from  A/D  converter  or  signal 
processor  doesn't  have  the  exactly  same  interval.  Delta  time(dt)  is  given  as 
follow, 

dt  =  (total  record  length)  /  (number  of  data  point-1)  (A-2) 


3)  mean(dp,ain) 

This  subroutine  calculates  the  mean  value  and  removes  the  mean 
value  of  the  signal  under  the  condition  of  the  input  variable  mvopt.  In  the 
case  the  signal  has  the  DC  component,  it  is  recommended  to  remove  the 
mean  value.  A  highpass  digital  filter  must  use  for  the  elimination  of  DC 
component. 


4)  filter(m,b,t,bk) 

This  subroutine  generates  a  nonrecursive  (finite  impulse  response  - 
FIR)  filter  weights  and  uses  the  method  devised  by  Potter,  Bickford  and 
Glaze.  Nonrecursive  highpass  filter  was  used  for  the  elimination  of 
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undesired  low  frequency  components.  The  variables  m,  b,  t,  and  bk  are  the 
number  of  the  weights,  the  cutoff  frequency  of  the  filter,  the  time  interval 
and  the  storage  array  of  the  filter  weights,  respectively. 


The  basic  design  was  to  use  a  symmetric  filter  of  the  form. 


m 


y(i)  =  E  ^k  s(i-k) 

(A-3) 

k=-m 

with 

b  k  =  bk 

(A-4) 

and 

,  sin  2;r/;kt 

(A-5) 

bk  = 

Tlk 


where  b^.  is  the  filter  weights,  y(i)  is  the  filtered  signal,  s(i)  is  the  original 
signal,  b  is  the  cutoff  frequency,  t  is  a  sampling  interval  and  mo  is  the  span 
of  the  filter;  2m+l  weights  are  employed  because  of  symmetry,  only  m+1 
need  be  generated.  The  weights  are  computed  over  the  range  -m  to  m. 

The  weights  are  multiplied  by  a  window  function.  Potter  discusses  a 
number  of  windows  in  the  referenced  work.  His  P310  window  was  found  to 
be  appropriate  for  filter  implementation.  It  takes  the  form, 


=  —  do  +  2  X  d  cos^ 

M’  [  p^.3  m  _ 

where 

Cl.  =  —  k  =  ±  n 

k  2 

=  1  otherwise 

do  =  1 

d.j  =  di  =  0.684988 
d.2  =  d.  =  0.202701 
d.3  =  d3  =  0.0177127 


(A-6) 


(A-7) 


and 


3 

M-  =  do  +  2X  d  =  2.8108034  (A-8) 

p=-3 

For  a  highpass  filter  with  pass  band  from  the  cutoff  frequency(B)  to  the 
maximum  frequency,  generate  a  low  pass  filter  on  the  range  0  -  B,  and  then 
subtract  the  central  weight  from  unity  and  change  the  signs  of  the 
remainder  of  the  weights. 

The  filtered  signal  in  the  main  program  is  generated  by  using  the 
calcultaed  filter  weights  as  follows. 

do  k=-m,m 

j=k 

if  (k.lt.O)  then 
j=k*(-l) 
endif 


(A-9) 

bb=-bk(j) 
if  (k.eq.O)  then 
bb=l.-bk(j) 
endif 

y(i)=y(i)+bb*s(i-k) 
end  do 


5)  hammg(dp,dt,pi,ain) 

This  subroutine  is  a  data  tapering  by  using  the  modified  Hamming 
window.  It  is  often  desirable  to  taper  a  random  time  series  at  each  end  to 
enhance  certain  characteristics  of  the  spectral  estimates.  Tapering  is 
multiplying  the  time  series  by  a  data  window  analogous  to  multiplying  the 
correlation  function  by  a  lag  window.  Thus  tapering  the  time  series  is 
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equivalent  to  applying  a  convolution  operation  to  the  raw  Fourier 
transform.  The  purpose  of  tapering  when  viewed  from  its  frequency 
domain  effect  is  to  suppress  large  side  lobes  in  effective  filter  obtained  with 
the  raw  transform.  When  looked  at  from  the  time  domain,  the  object  of 
tapering  is  to  round  off  potential  discontinuities  at  each  end  of  the  finite 
segment  of  the  time  history  being  analyzed. 

The  used  data  window'  in  this  program  is  a  modified  Hamming 
wdndow'  as  given  in  Eq.(A-lO). 

0.54 -0.46*  cos(10;rt/T)  0<t<T/10. 

W(t)={1.0  T/10<t<9T/10  (A-10) 

0.54  -  0.46  *  cos(  107i(T-t)/T)  9T/10  <  t  <  T 


6)  anal(dp,pi,ain,s) 


This  subroutine  converts  a  real  signal  to  an  analytic  one  by  using  the 
Hilbert  transform  given  as  follow, 


h(n) 


(A-11) 


2  sin“(;m/2) 
n  n 
0, 


n  T!:  0, 
n  =  0. 


The  Wigner  distribution  has  a  periodicity  of  N/2,  where  N  is  the 
number  of  data.  Hence,  even  when  the  sampling  of  signal  satisfies  the 
Nyquist  criteria,  there  are  still  aliasing  components  in  the  Wigner 
distribution  function.  If  we  sample  the  signal  with  Nyquist  rate,  its  power 
density  spectrum  does  not  overlap  with  its  own  components.  For  Wigner 
distribution  function,  there  are  additional  spectrum  components  causing 
the  aliasing  interference.  However,  a  simple  way  to  alleviate  this  problem 
for  practical  purpose  is  to  increase  the  seperation  of  the  spectrum  groups  by 
either  doubling  sampling  rate  or  by  interpolating  w'ith  additional  data 


points  according  to  sampling  theorem  before  the  transform  process. 
Another  approach  to  avoid  the  aliasing  is  to  use  only  the  positive  part  the 
signal's  frequency  components,  the  analytical  signal,  before  computing  the 
Wigner  distribution  function.  The  analj'tic  signal  approach  is  convenient 
since  it  can  be  easily  obtained  with  the  Hilbert  transform,  which  can  also 
use  the  efficient  FFT  algorithm. 

This  program  uses  the  later  case.  The  real  time  signal  converts  to 
complex  variable  by  the  Hilbert  transform,  that  is,  the  imaginary  part  is 
generated  by  the  convolution  of  the  impulse  response  h(n)  shown  as  follow', 

oc 

H(Sr(n)}  =  21  (A-13) 


where  Sr  is  the  original  input  signal. 

The  analytic  version  of  the  real  signal  is  made  up  of  the  real  signal 
plus  an  imaginary  part  composed  of  the  Hilbert  transform  of  the  real 
signal. 

S(l)  =  Sr(t)  +  j  H{Sr(t)}  (A-14) 


7)  wigner(dp,dt,pi,s,c,mm,nn,wdf,mt,nO 

This  subroutine  calculates  the  Wigner  distribution  function.  The 
Wigner  distribution  function  is  obtained  by  using  fast  Fourier  transform 
(FFT)  of  the  local  correlation  of  sigTsal.  At  first,  the  local  correlation  of  the 
signal  is  obtained  and  2N  point  :  FFT  is  carried  out.  This  program 
generates  the  results  wdfli,j)  of  the  Wigner  distribution  function  and 
pseudo  Wigner-Ville  distribution  function  on  the  given  time  series.  The 
Wigner  distribution  has  1/2  of  the  frequency  resolution  of  the  power 
spectrum  because  2N  points  FFT  and  the  argument  of  the  time  signal  and 
its  conjugate  contains  a  factor  of  1/2.  That  is,  df  =  l/(4*dp*dt). 


To  avoid  the  negative  values  and  the  elimination  of  the  interference 
by  the  cross  correlation  term,  this  program  uses  the  smoothing  technique 
with  Gaussian  window  function  as  following  Eq.(A-15). 

t'  (O' 

G(t,w)  =  - — ! - e  (A-15) 

2;rcT,(T^ 

The  smoothing  is  obtained  by  the  convolution  integration  the  Wigner 
distribution  and  a  Gaussian  window’  function  as  given  Eq(26).  To  perform 
the  convolution  on  the  sampled  wdf(i,j),  the  Gaussian  function  w'as  first 
truncated  so  that  it  spanned  the  range  ±  2(T,  and  ±2  0“^.  To  avoid  the 

negative  value,  Gaussian  window  size  must  have  a  larger  or  equal  integer 
value  mtand  nf  than  (dp/7r)l/2  Por  example,  in  the  case  of  dp=1024,  at  least 
Gaussian  window  size  must  be  18x18  or  larger.  We  recommend  to  select 
the  square  windows  size. 


8)  ffudp,pi,c) 

This  subroutine  is  a  program  about  the  Fast  Fourier  Transform  by  Jim 
Cooley's  method. 


APPENDIX  B. 


PROGRAM  LIST 


B-1.  PROGRAXI  LIST  ABOUT  THE  CALCULATION  OF  WDF 


PROGRAM  PWA'D 

♦  * 

*  PSEUDO  WIGNER-VILLE  DISTRIBUTION  FUNCTION  ♦ 

*  pc  VERSION  1.0  JAN.  1993  * 

*  bv  I'sing  highpass  digital  filter)  * 

*  * 
:(c!)c*j|t*  +  *  +  4:*>f:  +  *  +  +  +  +  *  +  +  +  *4;!i<*****  +  **!it**=(<******>f*!(<**+ **************** 


* 

★ 

* 

* 

♦ 

* 

♦ 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 

* 

* 

* 


This  is  the  program  about  the  pseudo  Wigner-Ville 
distribution(PWVD).  The  P\WD  is  a  three  dimensionaKtime, 
frequency,  amplitude)  representation  of  an  input 
signal  and  is  ideally  suited  for  portrajdng  transient 
phenomena. 

VARIABLES 

dp  =  the  number  of  the  sampled  data  point 
mm  =  output  parameter  for  time 
nn  =  output  parameter  for  frequency 
mt  =  smoothing  parameter  for  time  in  Gaussian  Fnc. 
nf  =  smoothing  parameter  for  frequency 
in  Gaussian  Fnc. 

inname  =  input  filename  (sampled  data) 
mvopt  =  option  with  respect  to  remo\ing  mean  value 
if  mvopt=l,  then  zero  mean 
if  mvopt=0,  then  no 


df  =  frequency  resolution 

dt  =  sampling  interval 


*  ARRAYS 

♦ 

*  tin(i)  =  sampled  time  data 

*  ain(i)  =  sampled  magnitude  data 

*  fain(i)  =  filtered  data 

*  wdRij)  =  reduced  array  of  PWATD 

*  bk(i)  -  filter  weights 

*  s(i)  =  analytic  signal 

c(i)  =  local  auto-correlation  or  the  results  of  FFT 

* 

* 

*  VARIABLE  DECLARATION  * 

*  * 

Jtc***************:***************^******  **♦♦*)♦;♦*♦**♦♦***♦******♦***** 

number  of  data  points 

*  Notes:  For  sample  sizes  greater  than  2048,  change 

*  the  parameter  np. 


integer  dp,dp2,mvopt,redopt,mm,nn,nf,mt 
parameter  (np=2048) 

real  pi,  tin(np), ain(np),wdf\256, 128),  dt, faint  np), 

1  bk(llOO) 

complex  s(np*2),c(np*2) 
character*25  inname 
character  as 
pi=atan(l.)*4. 

********!♦;*♦  +  **  +  *  +  *  +  ******  +  *****  +  +  *******♦  +  **>(<****♦*****♦***  + 

*  —  Print  description  of  program  — 

print*  '******************♦*+♦♦♦**♦***♦**********♦*' 
print*,'  Pseudo  Wigner-Ville  Distribution’ 


print* 


***************  *********' 


print* 


--  Set  input  parameters  — 

print*,’  Enter  name  of  signal  input  file' 

read(5,901)  inname 

print* 

print*, 'Number  of  the  sampled  data  point’ 
read(5,902)  dp 

dp2=dp*2 

do  100  i=l,dp2 
s(i)=cmplx(0.,0.) 

print*,'  Do  you  vsTsh  to  remove  the  mean  value?' 
print*,'  Enter  1  for  yes  or  0  for  no’ 
read(5.902)  mvopt 
print* 

print*.  Do  you  want  to  apply  a  highpass  digital' 
print*, 'filter  to  the  original  data  ?  (Y/'N)' 
read(*,903)  as 

if  (as.eq.'Y'.or.as.eq.'y')  then 

print *,'Ent«r  the  cutoff  frequency  of 
print*, 'the  digital  highpass  filter  (in  Hz)  ' 
readi*,*)  bw 


print*, 'Input  the  desired  reduction  size' 

print*,’  input  1  for  64  by  32  ' 

print*,'  input  2  for  128  by  64  ’ 

print*,'  input  3  for  128  by  128  ' 

print*,'  input  4  for  256  by  128  ’ 

read(5,902)  redopt 

print* 

if  (redopt.eq.l)  then 
mm=dp/32 
nn=dp2/64 

elseif  (redopt. eq. 2)  then 
mm=dp^64 
nn=dp2/128 

elseif  (redopt. eq, 3)  then 
mm=dp/128 
nn=dp2/128 
else 

mm=dp/128 

nn=dp2''256 

endif 

print*,  Input  the  desired  smoothing  window  size 
print* 

print*,'  input  the  smoothing  parameter  for  frequency' 
print*,'  in  Gaussian  function' 
read(5,902)  nf 

print*,'  input  the  smoothing  parameter  for  time' 
print*,'  in  Gaussian  function' 
read(5,902)  mt 


*  The  calculation  part  of  the  program  * 


Read  the  sampled  data  file 


* 

open(4,  file=inname,status='old’) 
call  indata(dp,tin,ain) 
print*, 'finished  subr.  indata' 
close(4) 

*  calculate  the  mean  time  interval 

call  dtcalc(dp,tin,dt) 
print*, 'finished  subr.  dtcalc' 

if  (mvopt.eq.  1 )  then 
call  meanidp.ain) 
print*,  'finished  subr.  mean' 
endif 

*  Signal  modifications 

* 

*  Application  of  highpass  digital  filter 

if  (as.eq.'V'.or.as.eq.'y')  then 
mo=dp/2 

*  calculate  the  filter  weighting 

call  filtenmo,bw,dt,bk) 

*  pass  the  highpass  filter 

do  160  i=l,dp 
160  fain(i)=0. 

do  200  i=l,dp 


60 


do  170  k=-mo,mo 

j=k 

if  (k.lt.O)  then 
j=kn-l) 
endif 

j=j+l 

ll=i-k 

if  (ll.lt.l)  then 
ll=ll+dp 

elseif  (Ihgt-dpl  then 
ll=ll-dp 
endif 
bb=-bk(j) 
if  (k.eq.O)  then 
bb=l-bk<j) 
endif 

fain(i)=fain(i  )+bb*ain(ll) 

170  continue 

200  continue 

do  210  i=:l,dp 
210  ainii  )=fain(i) 

print*,  'finished  digital  filtering' 
endif 

*  Window  application(modified  hamming  windows) 

call  hammg(dp,dt,pi,ain) 
print*,  'finished  subr.  hammg' 

*  Conversion  of  real  signal  to  an  analytic  one 
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call  anali dp.pi,ain,s) 


print*,  finished  subr.  anal' 


open(9.rile='r\vdrout',status='ne\v’) 
open!  10,rile='rswdf.out',status='ne\v‘) 

♦ 

*  Writing  the  WDF  to  output  file  * 

ttime=dp*dt 
df=l./(4.*dp*dt) 
fmax=2.*dp*df 
nx=dp2/nn 
ny= dp/mm 

v,Tite(9,904)  tint l),ttime,fmin,fmax 
write(9,*)  nx,ny 

*  Calculation  of  the  Wigner  distribution: 

print*,  Calculating  the  PW\T)' 
call  wigner(dp,dt,pi,s,c,mm,nn.wdf,mt,nf) 
print*. 'finished  wigner  ville  distribution  function' 
close(9) 

* 

*  Writing  of  reduced  &  smoothed  \WD  to  output  file 

♦ 

writet  10,904)  tint  1  ),ttime,fmin,fmax 
writetlO,*)  nx,ny 
do  500  i  =  l,dp/nim 
do  500  j=l,dp2/nn 

writetl0,905)  wdflj.i) 

500  continue 

*  Format  statements 
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s  s  s  s 

cn  CO  to 


901  format(a25) 
format(i6) 
format(al) 

format(el2.5,/,el2.5/,el2.5/,el2.5) 
format(2x,el2.5) 

close(lO) 

print*, 'terminate  the  excution  of  Wigner-Ville  distribution’ 

stop 
end 

*  * 

♦  SUBROUTINES  * 

♦  ♦ 
^i:^^i^^li:t:*^4i****^*****************^********************************* 

*******;f;***%,f;**+*********;(i  ******♦***♦*. t ♦♦♦♦♦♦♦**»♦*******★***♦ 

subroutine  indata(dp,tin,ain) 

integer  dp 
real  tin(*),ain(*) 

********simple  loop  to  read  in  time  &  amplitude************* 

do  100  j=l,dp 

read(4,*)  tin(j)  ,  ain(j) 

100  continue 
return 
end 
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subroutine  dtcalc(dp,tin,dt) 

♦ 

*  This  subroutine  calculates  the  delta  t  of  the  signal 

* 

* 

integer  dp 
real  tin(*),dt 

♦ 

dtsum  =  0.0 

do  100  i  =  1 ,  dp-1 

delt  =  tin(i+l)  -  tin(i) 
dtsum  =  dtsum  +  delt 

100  continue 

dt  =  dtsum  /  float(dp-l) 

return 

end 

subroutine  mean(dp,ain) 

*  This  subroutine  calculates  and  removes  the  mean  value 

*  of  the  signal. 

* 

integer  dp 
real  ain(*),meanv 
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asum  =  0.0 
do  100  i  =  1  ,  dp 

asum  =  asum  +  ain(i) 
100  continue 

meanv  =  asum  /  dp 
do  200  i  =  1  ,  dp 

ain(i)  =  ain(i)  -  meanv 
200  continue 

return 

end 


subroutine  rilter(mo,b,t,bk) 

♦ 

*  Routine  generates  FIR  filter  weights. 

*  Method  devised  by  Potter,  Bickford  and  Glaze. 

*  There  are  a  total  of  2M+1  weights... filter  generates  M+1. 

* 

*  —  variables  — 

*  t  =  the  sampling  interval  in  second. 

*  bw  =  cutoffi half-power  point)  of  the  filter  in  Hz; 

*  must  be  on  the  range  from  o  to  l/2t. 

* 

*  Results  are  stored  in  bk 
+ 

*  Note;  in  the  case  of  highpass  filter,  the  value  of 

*  weight  bo  must  use  1-bO  instead  of  bO. 

dimension  bk(*),d(3) 

data  d0/0.35577019/,d(l)/0.2436983/,d(2)/0.07211497/, 

*  d{3)/0.00630165/ 
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pi=atan(l.)*4. 

m=mo 

*  first  generate  plain  boxcar  weights 
fact=2.*b*t 

bk(l)=fact 
fact=fact*pi 
do  5  i=l,m 
fi=i 

5  bk(i+l)=sin(fact*fi)/(pi*fi) 

*  trapezoidal  weighting  at  end 
bk(m+l)=bk(m+l)/2. 

*  Now  apply  the  Potter  p3 10  window 
sumg=bk(l) 

do  15  i=l,m 
suni=d0 

fact=pi*fioat(i)/float(m) 
do  10  k=l,3 

10  sum=sum+2.*d(k)*cos(fact*float(k)) 

bk(i+l)=bk(i+l)*sum 
15  sumg=sumg+2.*bk(i+l) 

ml=m+l 
do  20  i=l,nil 
20  bk(i)=bk(i)/sumg 

return 
end 


subroutine  hammg(dp,dt,pi,ain) 

♦ 

*  This  subroutine  applies  a  modified  hamming  window 

*  to  the  signal  ain(t) 

* 

integer  dp 
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real  pi,ain( *),dt,mtime,del  1  ,de!2, const 

* 

nitime=(dp-l)*dt 

dell=0.1*nitime 

del2=0.9*mtime 

const=piydell 

do  100  j  =  1  ,  dp 

t  =  (j-l)*dt 

if  (t.le.dell)  then 

ain(j)  =  ain(j)  *  (.54-0.46*cos(consi*t)) 
elseif  ((t.ge.del2).and.(t.le.mtime))  then 

ain(j)=ain(j)*(.54-0.46*cos(const*(mtime-t))) 

endif 

100  continue 

return 

end 


subroutine  anal(dp,pi,ain,s) 

♦  ♦  +  +  +  **+*  +  ********)(c***Jt:  +  +  *  +  *  +  *  +  *****************  +  ********=)‘ 

♦ 

*  This  subroutine  converts  a  real  signal  to  an 

*  analytic  one  by  using  Hilbert  transform. 

* 

♦ 


s(i)  =  analytic  signal. 


integer  dp 

real  pi,ain(*),sum,sumb,val,sval 
complex  s(*) 

do  100  i=l,dp 
sum=0.0 

do  200  j=l,dp 
sumb=0.0 

if  (i-j.eq.O)  go  to  200 


n=i-j 

val=pi*rL^2. 

sval=sin(val) 

sumb=ain(j)*sval*sval/val 

200  sum=sum+sumb 

s(i)=cmplx(ain(i),sum) 

100  continue 

return 

end 


subroutine  wigner(dp,dt,pi,s,c,mm,nn,wdf,mt,nO 

+  **  +  **  +  ******  +  *****3f:*3ic  +  5f:*  +  =t:+  +  +  %**a»'  +  +  *  +  +  +  **  +  +  +  *  +  +  +  *  +  +  +  +  *+  + 
* 

*  This  subroutine  calculates  the  WDF  of  the  signal 

* 

* 
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integer  dp,dp2 

real  pi, dt,coef,wdfl256, 128), hg(-60:60, -60:60) 

complex  s(*),dum,c(*) 

dp2  =  dp*2 

coef  =  2.0  *  dt 

df=l./(4.*dp*dt) 

nl2=nf*2 

mt2=mt*2 

fl=float(mt) 

f2=float(n0 

*  Gaussian  function 

* 

val=l./((2.*pi)**2*fl*f^*dPdt) 
do  20  j=:-mt2,mt2 
qi=float(j) 
do  10  i=-n02,nf2 
q2=float(i) 

cf=-((ql*ql)/(2.*n*fl))-((q2*q2)/(2.*f2*f2)) 

hg(i,j)=val*exp(cf) 

10  continue 

20  continue 

* 

*  initialize  wdflij) 

♦ 

do  100  i=  1,256 
do  100j=l,128 
100  wdflij)=0. 

do  5000  j  =  1  ,  dp 

*  local  auto-correlation 

do  1000  i  =  1 ,  dp+1 
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1000 


1400 

1500 


if  (j.ge.i)  then 
dum  =  s(j-i+l) 
else 

dum  =  cmplx(0.,0.) 
endif 

c(i)  =  coef  *  (s(j+i-l)*conjg(dum)) 
if  (i.ne.l.and.i.ne.dp+1)  then 
c(dp2-i+2)  =  conjg(c(i)) 
endif 
continue 

call  fft(dp,pi,c) 


ik=mod((j-l),mm) 
if  (ik.eq.O)  then 

do  1500  i=l,dp2,nn 

write(9,1400)  real(c(i)) 
format(2x,el2.5) 
continue 
endif 
ml=0 
do  4000 
m=l,dp,mm 
ml=ml+l 
nl=0 

if  (abs(j-m).le.mt2)  then 
do  3000  n=l,dp2,nn 
nl=nl+l 

do  2500  kk=n-nf2,n+nf2 
kl=kk 

iflkk.lt. 1)  kl=kk+dp2 
iflkk.gt.dp2)  kl=:kk-dp2 

wdf\nl,ml)=wdflnl,ml)+real(c(kl))*hg(kk-nJ-m)*dPdt 
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2500 

continue 

3000 

continue 

endif 

4000 

continue 

5000 

continue 

return 

end 


subroutine  fTt(dp,pi,c) 

* 

*  This  subroutine  is  the  Fast  Fourier  Transform 

* 

integer  dp,dp2,val,coef,coefl 
real  pi 

complex  dum,c(*),dum3,dum2 

* 

dp2  =  dp*2 

const=float(dp2) 

val=alog(const)/alog(2.)+.l 

j=l 

do  40  i=l,dp2-l 
if  (i.ge.j)  go  to  10 
dum3=c(j) 
c(j)=c(i) 
c(i)=dum3 
10  k=dp 

20  if  (k.ge.j)  go  to  30 
j=j-k 
k=k/2 
go  to  20 
30  j=j+k 
40  continue 


do  70  n=l,val 
coef=2**n 
coefl=coeC'2 
dum2=cmplx(  l.,0.) 
theta=pL'float(coefl ) 
dum=cmplx(cos(  theta  ),-sin(theta)) 
do  60  j=l,coefl 

do  50  i=j,dp2,coef 
ii=i+coefl 
dum3=c(ii)*du 
m2  c(ii)=c(i)- 
dum3 

c(i  )=c(i)+dum3 


50 

continue 

dum2=dum2 

60 

continue 

70 

continue 

return 

end 
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PROGRAM  LIST  OF  3-D  PLOT 


PROGRAM  3DPLOT 

This  program  uses  the  graphic  package  CA-DISSPLA  to  plot 
the  results  of  WTF 

ttime  =  time  record  length 

tini  =  initial  time 

fmin  =  start  frequency 

fmax  =  stop  frequency 

nx  =  the  number  of  the  frequency  data  at  time  n  x  dt 

ny  =  the  number  of  the  time  data  at  frequency  m  x  df 

.  Declaring  variables . 

real  rvvdT32768) 
integer  nx,ny 
character*25  fname 

writef*,*)  'input  file  name  ?' 
read(*,20)  fname 
format(a25) 

open(15,file=fname,status='old') 

read(15,*)  tini 
readdS,*)  ttime 
readdS,*)  fmin 
readdS,*)  fmax 
readd5,*)  nx,ny 

write!*,*)  tini, ttime, fmin, fmax, nx,ny 


n=nx*ny 
do  100  i=l.n 

readUS,*)  nvdfli) 

100  continue 

close(15) 

smax=rwdftl) 
do  200  i=l,n 

if  (smax.lt.rwdfu))  then 
smax=r\vdf(i) 
endif 

200  continue 

writei*,*)  'smax  =',  smax 

write!*,*)  'input  maximum  z-axis  value' 

read!*,*)  fac 

call  pdev('ln03',  ieer) 

*  plotting 

call  hwshd 
call  swissm 

call  shdchr!90., 1,0. 002,1) 
call  height(0.2) 
call  physor(0.7,0.625) 
call  area2d!7.5,9.75) 

call  messag!'WIGNER-VILLE  DISTRIBUTION  $',100,1.1,8.2) 
call  blsur 

call  volm3d!8.,8.,9.) 

call  x3name!'Frequency  !Hz)  S'.  100) 

call  ySnameCTime  (sec)  $',100.) 
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call  z3name('Amplitude  $',100) 

call  vuangl(-60.,30.,30.) 

call  zaxangOO.) 

tstep=ttime/4. 

fstep=(fma':-fmin)/4. 

t.max=tini+ttime 

call  graf3d(rimn,fstep,fmax,tini,tstep,tmax,0.,'SCALE',fac) 

call  surmat(rwdf,l.nx,l,ny,l) 

call  end3gr(0) 
call  endpl(o) 
call  donepl 

stop 

end 
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B-3.  PROGRA>I  LIST  OF  CONTOUR  PLOT 


PROGRAM  CONTOUR 


*  This  program  uses  the  graphic  package  CA-DISSPLA  to  plot 

*  the  results  of  \VT)F 

* 


♦ 

ttime 

=  time  record  length 

♦ 

tini 

=  initial  time 

fmin 

=  start  frequency 

♦ 

fmax 

=  stop  frequency 

* 

nx 

=  the  number  of  the  frequency  data  at  time  n  x  dt 

* 

ny 

=  the  number  of  the  time  data  at  frequency  m  x  df 

♦ 


*  .  Declaring  variables . 

* 

real  r(64,128) 
common  work(9000' 
integer  nx.nv 
character*25  fname 

write(*,*)  'input  file  name  ?' 
reach*, 20)  fname 

20  format(a25) 

* 

open(15,file=fname,status=’old') 

reach  15,*)  tini 
reachl5,*)  ttime 
reach  15,*)  fmin 
reachl5,*)  fmax 
reach  15,*)  nx,ny 

write!*,*)  tini, ttime, fmin, fmax, nx,ny 


smax=0. 
do  100  i=l,ny 

do  100  j-  l,nx 
read(.15,*)  r<ij) 
if  (r(iJ).lt.O.)  Kij)=0. 
if  (smax.lt.r(ij))  then 
smax=r(i,j) 
endif 

100  continue 

close(15) 

*  Normalizing 

do  200  i=l,ny 
do  200  j=l,iix 

r<ij)=r(ij)/sniax 
200  continue 

write(*,*)  'smax  smax 

writer*,*)  'input  maximum  z-axis  value’ 

readC*,*)  fac 

call  pdevClnOS',  ieer) 

*  plotting 

call  hwshd 
call  swissm 

call  shdchr(90., 1,0. 002,1) 
call  height(0.2) 
call  page(8.5,ll.) 
call  physor(1.5,1.5) 
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call  area2d(6.,6.) 

call  headinCWIGNER-VILLE  DISTRIBUTION  CONTOURS’, 
*  100,1.1,1) 
call  blsur 

call  ynameCFrequc  -icy  (Hz)  $’,100) 

call  xnameC’Time  f'sec)  $’,100) 

call  yaxangCO.) 

tstep=ttime/5. 

fstep=(fmax-fmin)/5. 

tmax=tini+ttiine 

call  graf(tini,tstep,tmax,fimn,fstep,fmax) 

call  frame 

call  bcomonOOOO) 

scale=l./30. 

call  conmak(r,ny,nx, scale) 
contour  plot 

call  conlin(0, ’SOLID  ,  ’NOLABELS’,1,3) 
call  conIin(l.,'DASH',  ’NOLABELS’, 1,1) 
call  conang(50.) 
call  raspln(0.25) 

call  contur(2, ’LABELS’, ’DRAW’) 

call  endpKo) 
call  donepl 

stop 

end 
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B-4.  PROGRAM  LIST  OF  PLOT  OF  TIME  SERIES 


PROGRAM  TIMEPLOT 

* 

*  This  program  uses  the  gi-aphic  package  CA-DISSPLA  to  plot 

*  the  time  series. 

*  .  Declaring  variables  - . 

* 

dimension  x(4000),y(4000) 
character*60  title 
character*25  inname 

write(*,*)  'input  file  name  ?' 
read(*,1000)  inname 

1000  formate  a25) 
write(*,*)  'title  ?' 
read(*,1001)  title 

1001  format(a60) 

write(*,*)  'the  number  of  the  data  point’ 
read(*,*)  np 

write(*,*)  'maximum  scale  of  the  time  in  a  figure' 
read(*,*)  xmax 

write(*,*)  'maximum  scale  of  the  magnitude  in  the  figure' 
read(*,*)  yrnax 

* 

open(8,file=inname,status='old') 

do  100  i=l,np 

read(8,*)  x(i),y(i) 

100  continue 
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close(8) 


call  pdevClnOS',  leer) 

plotting 

call  hwshd 
call  swissm 

call  shdchr{90., 1,0. 002,1) 

call  height(0.2) 

call  page(8.5,ll.) 

call  physor(1.5,1.5) 

call  area2d(6.,6.) 

call  xnameCTime  (sec)  $’,100) 

call  ynameCFrequency  (Hz)  $',100) 

call  headin( title, 60, 1.1,1) 

call  thkfrm(O.Ol) 

Ck'U  yaxangOO.) 

call  graf(0.,xmax/4.,xmax,-ymax,0.5,ymax) 

call  grid(l,l) 

call  curv'e(x,y,np,0) 

call  endpl(o) 
call  donepl 

stop 

end 
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